home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / ddistribs.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  10KB  |  438 lines
  1. /* ddistributions - Basic discrete probability distributions           */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. #include "xlisp.h"
  8. #include "osdef.h"
  9. #ifdef ANSI
  10. #include "xlproto.h"
  11. #include "xlsproto.h"
  12. #else
  13. #include "xlfun.h"
  14. #include "xlsfun.h"
  15. #endif ANSI
  16.  
  17. #ifdef LATTICE /* name conflict JKL */
  18. #define rand randl
  19. #endif
  20.  
  21. /* forward declarations */
  22. #ifdef ANSI
  23. LVAL binomialcdf(void),poissoncdf(void),binomialpmf(void),poissonpmf(void),
  24.      binomialquant(void),poissonquant(void),binomialrand(void),
  25.      poissonrand(void),cdf(char),pmf(char),quant(char),findcompound(int),
  26.      rand(char);
  27. void getbinargs(int *,double *),getpoisarg(double *),
  28.      fixuparglist(LVAL);
  29. double binomial_cdf(int,int,double), poisson_cdf(int,double);
  30. int binomial_quant(double,int,double),poisson_quant(double,double),
  31.     poisson_rand(double),binomial_rand(int,double),findrlen(LVAL);
  32. #else
  33. LVAL binomialcdf(),poissoncdf(),binomialpmf(),poissonpmf(),
  34.      binomialquant(),poissonquant(),binomialrand(),
  35.      poissonrand(),cdf(),pmf(),quant(),rand(),findcompound();
  36. void getbinargs(),getpoisarg(),fixuparglist();
  37. double binomial_cdf(), poisson_cdf();
  38. int binomial_quant),poisson_quant(),poisson_rand(),binomial_rand(),findrlen();
  39. #endif ANSI
  40.  
  41. #ifndef PI
  42. # define PI 3.14159265358979323846
  43. #endif PI
  44. #ifndef nil
  45. # define nil 0L
  46. #endif
  47.  
  48. /* numerical distribution function */
  49. LOCAL LVAL binomialcdf()  { return(cdf('B')); }
  50. LOCAL LVAL poissoncdf()   { return(cdf('P')); }
  51.  
  52. /* recursive distribution functions */
  53. LVAL xsrbinomialcdf()
  54.     { return(recursive_subr_map_elements(binomialcdf, xsrbinomialcdf)); }
  55. LVAL xsrpoissoncdf()
  56.     { return(recursive_subr_map_elements(poissoncdf, xsrpoissoncdf)); }
  57.  
  58. /* numerical probability mass function */
  59. LOCAL LVAL binomialpmf()  { return(pmf('B')); }
  60. LOCAL LVAL poissonpmf()   { return(pmf('P')); }
  61.  
  62. /* recursive probability mass functions */
  63. LVAL xsrbinomialpmf()
  64.     { return(recursive_subr_map_elements(binomialpmf, xsrbinomialpmf)); }
  65. LVAL xsrpoissonpmf()
  66.     { return(recursive_subr_map_elements(poissonpmf, xsrpoissonpmf)); }
  67.  
  68. /* numerical quantile function */
  69. LOCAL LVAL binomialquant()  { return(quant('B')); }
  70. LOCAL LVAL poissonquant()   { return(quant('P')); }
  71.  
  72. /* recursive probability mass functions */
  73. LVAL xsrbinomialquant()
  74.     { return(recursive_subr_map_elements(binomialquant, xsrbinomialquant)); }
  75. LVAL xsrpoissonquant()
  76.     { return(recursive_subr_map_elements(poissonquant, xsrpoissonquant)); }
  77.  
  78. /* random number generating functions */
  79. LOCAL LVAL binomialrand()  { return(rand('B')); }
  80. LOCAL LVAL poissonrand()   { return(rand('P')); }
  81.  
  82. /* recursive probability mass functions */
  83. LVAL xsrbinomialrand()
  84.     { return(recursive_subr_map_elements(binomialrand, xsrbinomialrand)); }
  85. LVAL xsrpoissonrand()
  86.     { return(recursive_subr_map_elements(poissonrand, xsrpoissonrand)); }
  87.  
  88. /* argument readers */
  89. static void getbinargs(pn,pp)
  90.     int *pn;
  91.     double *pp;
  92. {
  93.   LVAL Ln, Lp;
  94.   int n;
  95.   double p;
  96.   
  97.   Ln = xlgafixnum(); n = getfixnum(Ln);
  98.   Lp = xlgetarg(); p = makedouble(Lp);
  99.   xllastarg();
  100.   if (n <= 0) xlerror("n is too small", Ln);
  101.   if (p < 0.0 || p > 1.0) xlerror("p not between 0 and 1", Lp);
  102.  
  103.   if (pn != nil) *pn = n;
  104.   if (pp != nil) *pp = p;
  105. }
  106.  
  107. static void getpoisarg(plam)
  108.     double *plam;
  109. {
  110.   LVAL Llam;
  111.   double lam;
  112.   
  113.   Llam = xlgetarg(); lam = makedouble(Llam);
  114.   xllastarg();
  115.   if (lam < 0.0) xlerror("lambda is too small", Llam);
  116.   
  117.   if (plam != nil) *plam = lam;
  118. }
  119.  
  120. /* Numerical Cdf's */
  121. LOCAL LVAL cdf(dist)
  122.     char dist;
  123. {
  124.   LVAL x;
  125.   double p, dx, lam, dp;
  126.   int ix, n;
  127.  
  128.   x = xlgetarg(); dx = makedouble(x); ix = floor(dx);
  129.  
  130.   switch (dist) {
  131.   case 'B': getbinargs(&n, &p);
  132.             dp = binomial_cdf(ix, n, p);
  133.             break;
  134.   case 'P': getpoisarg(&lam);
  135.             dp = poisson_cdf(ix, lam);
  136.             break;
  137.   default:  xlfail(" unknown distribution");
  138.   }
  139.   
  140.   return(cvflonum((FLOTYPE) dp));
  141. }
  142.  
  143. /* Numerical Pmf's */
  144. LOCAL LVAL pmf(dist)
  145.     char dist; /* chnaged JKL */
  146. {
  147.   LVAL x;
  148.   double p, dx, lam, dp;
  149.   int ix, n;
  150.  
  151.   x = xlgafixnum(); ix = getfixnum(x); dx = ix;
  152.  
  153.   switch (dist) {
  154.   case 'B': getbinargs(&n, &p);
  155.             if (p == 0.0) dp = (ix == 0) ? 1.0 : 0.0;
  156.             else if (p == 1.0) dp = (ix == n) ? 1.0 : 0.0;
  157.             else if (dx < 0.0 || dx > n) dp = 0.0;
  158.             else {
  159.               dp = exp(gamma(n + 1.0) - gamma(dx + 1.0) - gamma(n - dx + 1.0)
  160.                        + dx * log(p) + (n - dx) * log(1.0 - p));
  161.             }
  162.             break;
  163.   case 'P': getpoisarg(&lam);
  164.             if (lam == 0.0) dp = (ix == 0) ? 1.0 : 0.0;
  165.             else if (dx < 0.0) dp = 0.0;
  166.             else {
  167.               dp = exp(dx * log(lam) - lam - gamma(dx + 1.0));
  168.             }
  169.             break;
  170.   default:  xlfail(" unknown distribution");
  171.   }
  172.   
  173.   return(cvflonum((FLOTYPE) dp));
  174. }
  175.  
  176. /* Numerical Quantiles */
  177. LOCAL LVAL quant(dist)
  178.     char dist; /* chnaged JKL */
  179. {
  180.   LVAL x;
  181.   double p, dx, lam;
  182.   int n, k;
  183.  
  184.   x = xlgetarg(); dx = makedouble(x);
  185.   if (dx <= 0.0 || dx >= 1.0) xlerror("probability not between 0 and 1", x);
  186.   
  187.   switch (dist) {
  188.   case 'B': getbinargs(&n, &p);
  189.             if (p == 0.0) k = 0;
  190.             else if (p == 1.0) k = n;
  191.             else k = binomial_quant(dx, n, p);
  192.             break;
  193.   case 'P': getpoisarg(&lam);
  194.             if (lam == 0.0) k = 0;
  195.             else k = poisson_quant(dx, lam);
  196.             break;
  197.   default:  xlfail(" unknown distribution");
  198.   }
  199.   
  200.   return(cvfixnum((FIXTYPE) k));
  201. }
  202.  
  203. /* Random Number Generators */
  204. LOCAL LVAL rand(dist)
  205.      char dist; /* chnaged JKL */
  206. {
  207.   LVAL M, sample, variate;
  208.   double p, lam;
  209.   int m, n, i;
  210.  
  211.   M = xlgafixnum(); m = getfixnum(M);
  212.   if (m <= 0) xlerror("non positive number of variates", M);
  213.   
  214.   xlstkcheck(2);
  215.   xlsave(sample);
  216.   xlsave(variate);
  217.   
  218. /*  sample = NIL; not needed JKL */
  219.   
  220.   switch (dist) {
  221.   case 'B': getbinargs(&n, &p);
  222.             for (i = 0; i < m; i++) {
  223.               variate = cvfixnum((FIXTYPE) binomial_rand(n, p));
  224.               sample = cons(variate, sample);
  225.             }
  226.             break;
  227.   case 'P': getpoisarg(&lam);
  228.             for (i = 0; i < m; i++) {
  229.               variate = cvfixnum((FIXTYPE) poisson_rand(lam));
  230.               sample = cons(variate, sample);
  231.             }
  232.             break;
  233.   default:  xlfail(" unknown distribution");
  234.   }
  235.  
  236.   xlpopn(2);
  237.   
  238.   return(sample);
  239. }
  240.  
  241. LOCAL double binomial_cdf(k, n, p)
  242.     int k, n;
  243.     double p;
  244. {
  245.   double da, db, dp;
  246.   int ia, ib;
  247.   
  248.   if (k < 0) dp = 0.0;
  249.   else if (k >= n) dp = 1.0;
  250.   else if (p == 0.0) dp = (k < 0) ? 0.0 : 1.0;
  251.   else if (p == 1.0) dp = (k < n) ? 0.0 : 1.0;
  252.   else {
  253.   da = k + 1;
  254.   db = n - k;
  255.   ia = floor(da); ib = floor(db);
  256.   betabase(&p, &da, &db, &ia, &ib, &dp);
  257.    dp = 1.0 - dp;
  258.   }
  259.   return(dp);
  260. }
  261.  
  262. LOCAL double poisson_cdf(k, L)
  263.     int k;
  264.     double L;
  265. {
  266.   double dp, dx;
  267.   
  268.   if (k < 0) dp = 0.0;
  269.   else if (L == 0.0) dp = (k < 0) ? 0.0 : 1.0;
  270.   else {
  271.     dx = k + 1.0;
  272.     gammabase(&L, &dx, &dp);
  273.     dp = 1.0 - dp;
  274.   }
  275.   return(dp);
  276. }
  277.  
  278. LOCAL int binomial_quant(x, n, p)
  279.     double x, p;
  280.     int n;
  281. {
  282.   int k, k1, k2, del, ia;
  283.   double m, s, p1, p2, pk;
  284.   
  285.   m = n * p;
  286.   s = sqrt(n * p * (1 - p));
  287.   del = max(1, (int) (0.2 * s));
  288.   
  289.   k = m + s * ppnd(x, &ia);
  290.   k1 = k; k2 = k;
  291.   
  292.   do {
  293.     k1 = k1 - del; k1 = max(0, k1);
  294.     p1 = binomial_cdf(k1, n, p);
  295.   } while (k1 > 0 && p1 > p);
  296.   if (k1 == 0 && p1 >= x) return(k1);
  297.   
  298.   do {
  299.     k2 = k2 + del; k2 = min(n, k2);
  300.     p2 = binomial_cdf(k2, n, p);
  301.   } while (k2 < n && p2 < x);
  302.   if (k2 == n && p2 <= x) return(k2);
  303.   
  304.   while (k2 - k1 > 1) {
  305.     k = (k1 + k2) / 2;
  306.     pk = binomial_cdf(k, n, p);
  307.     if (pk < x) { k1 = k;/* p1 = pk; not used JKL */ }
  308.     else { k2 = k;/* p2 = pk; not used JKL */ }
  309.   }
  310.   return(k2);
  311. }
  312.  
  313. LOCAL int poisson_quant(x, L)
  314.     double x, L;
  315. {
  316.   int k, k1, k2, del, ia;
  317.   double m, s, p1, p2, pk;
  318.   
  319.   m = L;
  320.   s = sqrt(L);
  321.   del = max(1, (int) (0.2 * s));
  322.   
  323.   k = m + s * ppnd(x, &ia);
  324.   k1 = k; k2 = k;
  325.   
  326.   do {
  327.     k1 = k1 - del; k1 = max(0, k1);
  328.     p1 = poisson_cdf(k1, L);
  329.   } while (k1 > 0 && p1 > x);
  330.   if (k1 == 0 && p1 >= x) return(k1);
  331.   
  332.   do {
  333.     k2 = k2 + del;
  334.     p2 = poisson_cdf(k2, L);
  335.   } while (p2 < x);
  336.   
  337.   while (k2 - k1 > 1) {
  338.     k = (k1 + k2) / 2;
  339.     pk = poisson_cdf(k, L);
  340.     if (pk < x) { k1 = k;/* p1 = pk; not used JKL */ }
  341.     else { k2 = k; /* p2 = pk; not used JKL */ }
  342.   }
  343.   return(k2);
  344. }
  345.  
  346. /* poisson random generator from Numerical Recipes */
  347. LOCAL int poisson_rand(xm)
  348.     double xm;
  349. {
  350.   static double sqrt2xm, logxm, expxm, g, oldxm = -1.0;
  351.   double t, y;
  352.   int k;
  353.   
  354.   if (xm < 12.0) {
  355.     if (xm != oldxm) { expxm = exp(-xm); oldxm = xm; }
  356.     k = -1;
  357.     t = 1.0;
  358.     do {
  359.       k++;
  360.       t *= uni();
  361.     } while (t > expxm);
  362.   }
  363.   else {
  364.     if (xm != oldxm) {
  365.       oldxm = xm;
  366.       sqrt2xm = sqrt(2.0 * xm);
  367.       logxm = log(xm);
  368.       g = xm * logxm - gamma(xm + 1.0);
  369.     }
  370.     do {
  371.       do {
  372.         y = tan(PI * uni());
  373.         k = floor(sqrt2xm * y + xm);
  374.       } while (k < 0);
  375.       t = 0.9 * (1.0 + y * y) * exp(k * logxm - gamma(k + 1.0) - g);
  376.     } while (uni() > t);
  377.   }
  378.   return (k);
  379. }
  380.  
  381. /* binomial random generator from Numerical Recipes */
  382. LOCAL int binomial_rand(n, pp)
  383.     int n;
  384.     double pp;
  385. {
  386.   int j, k;
  387.   static int nold = -1;
  388.   double am, em, g, p, sq, t, y;
  389.   static double pold = -1.0, pc, plog, pclog, en, oldg;
  390.   
  391.   p = (pp <= 0.5) ? pp : 1.0 - pp;
  392.   
  393.   am = n * p;
  394.   if (p == 0.0) k = 0;
  395.   else if (p == 1.0) k = n;
  396.   else if (n < 50) {
  397.     k = 0;
  398.     for (j = 0; j < n; j++) if (uni() < p) k++;
  399.   }
  400.   else if (am < 1.0) {
  401.     g = exp(-am);
  402.     t = 1.0;
  403.     k = -1;
  404.     do {
  405.       k++;
  406.       t *= uni();
  407.     } while (t > g);
  408.     if (k > n) k = n;
  409.   }
  410.   else {
  411.     if (n != nold) {
  412.       en = n;
  413.       oldg = gamma(en + 1.0);
  414.       nold = n;
  415.     }
  416.     if (p != pold) {
  417.       pc = 1.0 - p;
  418.       plog = log(p);
  419.       pclog = log(pc);
  420.       pold = p;
  421.     }
  422.     sq = sqrt(2.0 * am * pc);
  423.     do {
  424.       do {
  425.         y = tan(PI * uni());
  426.         em = sq * y + am;
  427.       } while (em < 0.0 || em >= en + 1.0);
  428.       em = floor(em);
  429.       t = 1.2 * sq * (1.0 + y * y)
  430.         * exp(oldg - gamma(em + 1.0) - gamma(en - em + 1.0)
  431.               + em * plog + (en - em) * pclog);
  432.     } while (uni() > t);
  433.     k = em;
  434.   }
  435.   if (p != pp) k = n - k;
  436.   return(k);
  437. }
  438.